HW3_ExploratoryDataAnalysis_MXBJ

Dec 7 2018

MXBJ

# Installation of flowcatchR and calling necessary libraries
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("flowcatchR", version = "3.8")
## Bioconductor version 3.8 (BiocManager 1.30.4), R 3.5.1 (2018-07-02)
## Installing package(s) 'flowcatchR'
## 
## The downloaded binary packages are in
##  /var/folders/py/dxt709rd1nv01sr4bvsvbxrwvsqq8w/T//Rtmpkgtfee/downloaded_packages
## Update old packages: 'jsonlite', 'markdown'
library("flowcatchR")
## Loading required package: EBImage
library("ggplot2")
# Load image
# pre-drug image
fullData_predrug_long <- read.Frames(image.files="/Users/mbj268/Documents/PDA/Project/predrug-sp1-cropped.tif", nframes=1) 
## Creating a new object of class Frames...
## Created a Frames object of 1 frames.
inspect.Frames(fullData_predrug_long, nframes=36, display.method="raster")

# post-drug image
fullData_postdrug_long <- read.Frames(image.files="/Users/mbj268/Documents/PDA/Project/postdrug-sp3-cropped-long.tif", nframes=1)
## Creating a new object of class Frames...
## Created a Frames object of 1 frames.
inspect.Frames(fullData_postdrug_long, nframes=36, display.method="raster")

# Obtain images of the red channel
# pre-drug image
redchannel_predrug_long <- channel.Frames(fullData_predrug_long, mode="red")
inspect.Frames(redchannel_predrug_long, nframes=36, display.method="raster")

# post-drug image
redchannel_postdrug_long <- channel.Frames(fullData_postdrug_long, mode="red")
inspect.Frames(redchannel_postdrug_long, nframes=36, display.method="raster")

# Increasing the signal-to-noise 
# pre-drug image
preprocessed_predrug_long <- preprocess.Frames(redchannel_predrug_long,
                                    brush.size=2, brush.shape="disc",
                                    at.offset=0.2, at.wwidth=7, at.wheight=7,
                                    kern.size=3, kern.shape="disc",
                                    ws.tolerance=1, ws.radius=1)
## Warning in makeBrush(brush.size, brush.shape, step = FALSE): 'size' was
## rounded to the next odd number: 3
inspect.Frames(preprocessed_predrug_long, nframes=36, display.method="raster")

# post-drug image
preprocessed_postdrug_long <- preprocess.Frames(redchannel_postdrug_long,
                                    brush.size=3, brush.shape="disc",
                                    at.offset=0.06, at.wwidth=7, at.wheight=7,
                                    kern.size=3, kern.shape="disc",
                                    ws.tolerance=1, ws.radius=1)

inspect.Frames(preprocessed_postdrug_long, nframes=36, display.method="raster")

# Extract particles of interest
# pre-drug image
sampleparticles_predrug_long <- particles(redchannel_predrug_long, preprocessed_predrug_long)
## Computing features in parallel...
## Done!
# post-drug image
sampleparticles_postdrug_long <- particles(redchannel_postdrug_long, preprocessed_postdrug_long)
## Computing features in parallel...
## Done!
# Visualize selected particles 
# pre-drug image
painted_predrug_long <- add.contours(raw.frames=fullData_predrug_long,
                                 binary.frames=preprocessed_predrug_long,
                                 mode="particles")
inspect.Frames(painted_predrug_long, nframes=36, display.method="raster")

# post-drug image
painted_postdrug_long <- add.contours(raw.frames=fullData_postdrug_long,
                                 binary.frames=preprocessed_postdrug_long,
                                 mode="particles")
inspect.Frames(painted_postdrug_long, nframes=36, display.method="raster")

# Extract area, x coordinate and y coordinate of particles detected
# Organize image data and initialize variables for the loop
parts<-list()
parts[[1]]<-sampleparticles_predrug_long
parts[[2]]<-sampleparticles_postdrug_long
filtered_area<-list()
filtered_radius<-list()
filtered_eccentricity<-list()

# 1st for loop: for all image sequences
 for (k in 1:2) {
# Initialize variables 
  radius <- matrix(data=NA, nrow=length(parts[[k]]), ncol=50)
  eccentricity <- matrix(data=NA, nrow=length(parts[[k]]), ncol=50)
  xcoord <- matrix(data=NA, nrow=length(parts[[k]]), ncol=50)
  ycoord <- matrix(data=NA, nrow=length(parts[[k]]), ncol=50)
  area <- matrix(data=NA, nrow=length(parts[[k]]), ncol=50)
  
    # 2nd for loop: for all frames
  for(i in 1:length(parts[[k]])) {  
    if(length(parts[[k]][[i]][[1]])>0){
    radius[i, 1:length(parts[[k]][[i]][[8]])] <- parts[[k]][[i]][[8]] 
    eccentricity[i, 1:length(parts[[k]][[i]][[4]])] <- parts[[k]][[i]][[4]] 
    xcoord[i, 1:length(parts[[k]][[i]][[1]])] <- parts[[k]][[i]][[1]] 
    ycoord[i,1:length(parts[[k]][[i]][[2]])] <- parts[[k]][[i]][[2]]
    area[i,1:length(parts[[k]][[i]][[6]])] <- parts[[k]][[i]][[6]] 
    }
    else{
    radius[i] <- 0 
    eccentricity[i] <-0
    xcoord[i] <- 0
    ycoord[i] <- 0
    area[i] <- 0
    }}
 
 # Keep the max for each frame 
  r <- matrix(data=NA, nrow=length(parts[[k]]), ncol=1)
  e <- matrix(data=NA, nrow=length(parts[[k]]), ncol=1)  
  x <- matrix(data=NA, nrow=length(parts[[k]]), ncol=1)
  y <- matrix(data=NA, nrow=length(parts[[k]]), ncol=1)
  A <- matrix(data=NA, nrow=length(parts[[k]]), ncol=1)
  
  for (j in 1:length(area)) { 
  r[j] <- max(radius[j])
  e[j] <- max(eccentricity[j])
  A[j]<- max(area[j])
  x[j]<- max(xcoord[j])
  y[j]<- max(ycoord[j])
  }

  # Filter frames by area: keep particles that fall within a threshold
  area_indx<-which(A>10)
  area_r1<-r[area_indx]
  area_e1<-e[area_indx]
  area_A1<-A[area_indx]
  area_x1<-x[area_indx]
  area_y1<-y[area_indx]
  
  filtered_area[[k]]<-data.frame(area=area_A1, x=area_x1, y=area_y1, ex=area_e1, r=area_r1 )
  
    # Filter frames by radius: keep particles that fall within a threshold
  radius_indx<-which(r>1)
  radius_r1<-r[radius_indx]
  radius_e1<-e[radius_indx]
  radius_A1<-A[radius_indx]
  radius_x1<-x[radius_indx]
  radius_y1<-y[radius_indx]
  
  filtered_radius[[k]]<-data.frame(area=radius_A1, x=radius_x1, y=radius_y1, ex=radius_e1, r=radius_r1 )
  
    # Filter frames by eccentricity: keep particles that fall within a threshold
  eccentricity_indx<-which(e>0.7)
  eccentricity_r1<-r[eccentricity_indx]
  eccentricity_e1<-e[eccentricity_indx]
  eccentricity_A1<-A[eccentricity_indx]
  eccentricity_x1<-x[eccentricity_indx]
  eccentricity_y1<-y[eccentricity_indx]
  
  filtered_eccentricity[[k]]<-data.frame(area=eccentricity_A1, x=eccentricity_x1, y=eccentricity_y1, ex=eccentricity_e1, r=eccentricity_r1 )
  

 } 
# View histogram of x position to eliminate outliers  

pre_x_area<-filtered_area[[1]][[2]]
pre_y_area<-filtered_area[[1]][[3]]
pre_x_radius<-filtered_radius[[1]][[2]]
pre_y_radius<-filtered_radius[[1]][[3]]
pre_x_eccentricity<-filtered_eccentricity[[1]][[2]]
pre_y_eccentricity<-filtered_eccentricity[[1]][[3]]

post_x_area<-filtered_area[[2]][[2]]
post_y_area<-filtered_area[[2]][[3]]
post_x_radius<-filtered_radius[[2]][[2]]
post_y_radius<-filtered_radius[[2]][[3]]
post_x_eccentricity<-filtered_eccentricity[[2]][[2]]
post_y_eccentricity<-filtered_eccentricity[[2]][[3]]

# View histogram
# Pre-drug, x position 
# Filter x position by area
ggplot(data = data.frame(pre_x_area)) +
    geom_histogram(mapping = aes(x = pre_x_area))+ labs(x="x position")+ggtitle("Pre-drug filtered by area")+theme(plot.title = element_text(hjust = 0.5))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

 # Filter x position by radius
ggplot(data = data.frame(pre_x_radius)) +
    geom_histogram(mapping = aes(x = pre_x_radius))+ labs(x="x position")+ggtitle("Pre-drug filtered by radius")+theme(plot.title = element_text(hjust = 0.5))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

 # Filter x position by eccentricity
ggplot(data = data.frame(pre_x_eccentricity)) +
    geom_histogram(mapping = aes(x = pre_x_eccentricity))+ labs(x="x position")+ggtitle("Pre-drug filtered by eccentricity")+theme(plot.title = element_text(hjust = 0.5))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

# Post-drug, x position 
 # Filter x position by area
ggplot(data = data.frame(post_x_area)) +
    geom_histogram(mapping = aes(x = post_x_area))+ labs(x="x position")+ggtitle("Post-drug (5 min) filtered by area")+theme(plot.title = element_text(hjust = 0.5))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

 # Filter x position by radius
ggplot(data = data.frame(post_x_radius)) +
    geom_histogram(mapping = aes(x = post_x_radius))+ labs(x="x position")+ggtitle("Post-drug (5 min) filtered by radius")+theme(plot.title = element_text(hjust = 0.5))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

 # Filter x position by eccentricity
ggplot(data = data.frame(post_x_eccentricity)) +
    geom_histogram(mapping = aes(x = post_x_eccentricity))+ labs(x="x position")+ggtitle("Post-drug (5 min) filtered by eccentricity")+theme(plot.title = element_text(hjust = 0.5))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

  # xy positions 

# Pre-drug
# Filter xy position by area
xy_pre_area<-data.frame(x=pre_x_area,y=pre_y_area)  
ggplot(data = data.frame(xy_pre_area)) + geom_point(mapping = aes(x = x, y = y))+ggtitle("Pre-drug positions filtered by area ")+theme(plot.title = element_text(hjust = 0.5)) 

 # Filter x position by radius
xy_pre_radius<-data.frame(x=pre_x_radius,y=pre_y_radius)  
ggplot(data = data.frame(xy_pre_radius)) + geom_point(mapping = aes(x = x, y = y))+ggtitle("Pre-drug positions filtered by radius ")+theme(plot.title = element_text(hjust = 0.5)) 

 # Filter x position by eccentricity
xy_pre_eccentricity<-data.frame(x=pre_x_eccentricity,y=pre_y_eccentricity)  
ggplot(data = data.frame(xy_pre_eccentricity)) + geom_point(mapping = aes(x = x, y = y))+ggtitle("Pre-drug positions filtered by eccentricity ")+theme(plot.title = element_text(hjust = 0.5)) 

# Post-drug
# Filter xy position by area
xy_post_area<-data.frame(x=pre_x_area,y=pre_y_area)  
ggplot(data = data.frame(xy_pre_area)) + geom_point(mapping = aes(x = x, y = y))+ggtitle("Post drug (5 min) positions filtered by area ")+theme(plot.title = element_text(hjust = 0.5)) 

 # Filter x position by radius
xy_post_radius<-data.frame(x=pre_x_radius,y=pre_y_radius)  
ggplot(data = data.frame(xy_pre_radius)) + geom_point(mapping = aes(x = x, y = y))+ggtitle("Post drug (5 min) positions filtered by radius ")+theme(plot.title = element_text(hjust = 0.5)) 

 # Filter x position by eccentricity
xy_post_eccentricity<-data.frame(x=pre_x_eccentricity,y=pre_y_eccentricity)  
ggplot(data = data.frame(xy_pre_area)) + geom_point(mapping = aes(x = x, y = y))+ggtitle("Post drug (5 min) positions filtered by eccentricity ")+theme(plot.title = element_text(hjust = 0.5))